
\documentclass{article}
\usepackage{verbatim}

\begin{document}
\section{Introduction}
The simulation of incompressible flows requires the solution of elliptic
problems in order to keep the flow divergence free as time
progresses. The elliptic problem can be written formally as
\begin{equation}
{\bf u}={\bf u}^*-{\rm grad}\phi,\,\,{\rm div}{\bf u}=0,
\end{equation}
 with ${\bf u}\cdot{\bf n}$ assigned on the boundary of the domain. Here ${\bf
 u}^*$ is given, and the problem consists in solving for $\phi$ in order to
 obtain ${\bf u}$. If the equation above is dicretized using a finite volume
 approach, it is natural to define the D components of the velocity on a
 staggered grid, with the component $u_i$ defined on the corresponding edge of
 the computational cell, and the potential $\phi$ defined on the cell
 centers. Assuming curvilinear coordinates for generality, such that the
 computational domain is a square box, the problem can be
 written as 
\begin{equation}
\partial_i\sqrt(g)g^{ij}\partial_j\phi=\partial_i\sqrt(g)u^{*i},
\label{eq:master}
\end{equation}
with $g^{ij}\partial_j\phi$ assigned on the boundaries, and $u^{*i}$ denoting
the contravariant components of the vector ${\bf u}^*$.

\section{Discretization issues}
 
The goal is to solve the eq.~\ref{eq:master} so that  ${\rm div}{\bf u}=0$ to
machine precision. Thus, the gradient, divergence and laplacian ($\equiv$
div(grad)) operators have to be discretized consistently. For Cartesian
geometries, the latter yields the standard second-order accurate 7-point 
(5-point) stencil in 3D (2D). The stencil becomes  quickly more complex
 for curvlinear coordinates, where the general second-order 3D stencil uses 19
 points. Also, the extrapolation of the potential to the ghostpoints is not
 trivial when the metric is not diagonal at the boundaries. Existing general
 purpose elliptic solvers such as MUDPACK do not handle problems like
 eq.~\ref{eq:master} well because the solution they provide only gives ${\rm
   div}{\bf u}=O(\Delta^P)$, P being the order of the discretization. Hence
 the need of a dedicated kit.

\section{Solution strategy}

\begin{itemize}
\item Generate a set of discretized gradient, divergence and Laplacian
  operators expressed as sparse matrixes written in compact storage.
\item Leverage existing sparse matrix solvers libraries and optmized BLAS. 
The kit uses PARDISO as a
  direct solver, and the set of preconditioned iterative solvers implemented
  in SPARSKIT2. Matrix vector multipications are handled using the efficient
  MKL implementation of BLAS routines, while sparse-matrix/sparse-matrix
  operations are based on SPARSKIT2 (waiting for INTEL to catch up).
\item Use object oriented approach as much as possible, hiding the datails
  from the end-user. 
\item Even though the standard projection problem requires Neumann boundary 
  conditions, Dirichlet are implemented as well. Thus, the same package can be
  used to solve parabolic problems as well.  
\end{itemize}

\section{Implementation}

The kit contains 4 modules packaged in a library plus .mod files. Refers to
the examples for details.
\subsection{ElliDef}
This module contains definition of the structures used by the kit. The
interface to the user is the TYPE(Elli) structure. The user has to set ndim (2
or 3) the dimensionality of the problem, mX (the
number of cell centers along each direction), dX (the grid spacing) and
mthbc(1:6)  which specifies the type of BC along the boundary i
(convention is i=1 is left, i=2 is right,  i=3 front (or bottom in 2D), i=4
back (or top in 2D), i=5 bottom (not used in 2D) and i=6 top (not used in
2D)). mthbc(i)=0 means neumann, mthbc(i)=InternalBC means Dirichlet.
The pointer bc points to a structure containing the boundary conditions to the
problem. The user sets bc(i)\%p. If i is Neumann (NM), then bc(i)\%p contains $u^*-u$ on
that boundary. IF the boundary is Dirichlet (DR), then bc(i)\%p contains the value
of $\phi$ on the boundary. Note that for NM conditions bc(i)\%p(1:nx,1:ny)
where nx and ny are the dimensions of the boundary. If the boundary is DR,
then has to be assigned on bc(i)\%p(0:nx+1,0:ny+1).

div contains the divergence, and p the solution. pold is used to provide an
initial guess if the iterative solver is used. The dimensions are
div(1:mX(1),1:mX(2),1:mX(3)), pold same as div and p(0:mX(1),0:mX(2),0:mX(3)).

qe(1:mX(1)+1,1:mX(2)+1,1:mX(3)+1,2 or 3) contains the edge velocities, which
have to be set by the user after ellikit is initialized (and memory made
available).  

The \begin{verbatim}TYPE(CRS_P)\end{verbatim} 
pointer G is an array of pointers to the div, grad and
Laplacian operators. In case the user wants to use them (or need them) they
are defined as follows: G(i)\%crsptr points to gradient in i direction (1:3),
div in i direction (4:6), Laplacian \begin{verbatim}(CSR_LAPL)\end{verbatim} 
and preconditioner
\begin{verbatim} (CSR_PRECOND). \end{verbatim}
 In view of the symmetry, the last two are contain only the
upper triangular part. 
 
Note that all the automatic storage
 required by the Elli structure is allocate and deallocate by the kit when it
 is initialized and finalized. No user intervention is required.

\subsection{Basic-matrix}

This module contains subroutines that generate the compact storage
representation of the boulding blocks of the operators. They are basically
left and right differences (and averages), and a bunch of indentity/extension
perators. 

\subsection{matrix-operators}

This module contains two subsets of routines. There are routines to perform
matrix/matrix, matrix/vector operations and related tasks on compact storage
sparse matrices and routines to generate the gradient operator and handle the
boundary conditions. The first set of routines implements an object-oriented
wrapper to the SPARSKIT2 and MKL routines that do the actual dirty job. For
example, if A, B and C are pointers of TYPE(CRS) and B and C
contain the representation of
two matrixes, than\begin{verbatim} A=>ADB(B,C)\end{verbatim}
 returns a pointer to their product (in sparse
storage). The gradient operator is defined to act on p (defined on cell
centers including ghost layer) and returns the gradient on the corresponding
cell edges. The ``divergence'' operator acting along direction i 
acts on the cell edges
and returns the value in the cell centers. Thus, the application of div(grad)
in direction i reduces the size of
2. \begin{verbatim}EXTRAPOLATE_BC\end{verbatim} 
provides a convenienent
way of extrapolating the pressure on the ghost points based on the pressure in
the interior, when bc conditions are NM). 

\subsection{elliptic-solvers}

This module contains the routines that a typical user will use to set up and
solve the problem. The workflow is as follows:
\begin{enumerate}
\item Enter the relevant information (problem size,dimension,...) in the Elli
  structure.
\item call \begin{verbatim}err=CREATE_OPERATORS(Elli).\end{verbatim}
 This function allocates spaces and
  defines all the operators.
\item For safety, check that the logical function \begin{verbatim}check(Elli,CHECK_BC)\end{verbatim}
  returns true.
\item For safety, zeros the boundary condition vectors.
\item Load the the edge velocities to be projected in qe.
\item initialize the Poisson solver using
  \begin{verbatim}POISSON_SOLVERS(Elli,INIT,solver,OP),\end{verbatim}
 where solver denotes the type
  of solver desired. The list is \begin{verbatim}PARDISO_SOLVER,IT_DBCG,IT_BCG,IT_CG,&
     IT_CGNR,IT_FOM,IT_TFQMR,IT_BCGSTAB,IT_DQGMRES,IT_GMRES
,IT_FGMRES.\end{verbatim}
 There are two type of solvers. The direct solver
  \begin{verbatim}(PARDISO_SOLVER)\end{verbatim}
 uses
  pardiso to generate the LU decomposition of the laplacian, to be used later
  for the solution. All the others are preconditioned iterative solvers. See
  SPARSKIT2 for details. OP is an optional argument. If the direct solver is
  chosen, than leave it out. If an iterative solver is chosen, the
  presence/absence of it determines the type of preconditioner used. IF
  \begin{verbatim}OP=CSR_PRECOND\end{verbatim}
 is used, than the kit uses as a preconditioner the
  tridiagonal solver extracted from the laplacian picking the offdiagonal with
  the largest values (assumed coming from the derivative along the last
  direction). IF an iterative solver is chosen, and OP is not specified, then
  an incomplete LU will be used (see SPARSKIT2 ilut routine). The parameters
  for the ilut are set in setparm. The user has the
  flexibility of providing a different preconditioner in compact storage,
  stored in \begin{verbatim}Elli\%G(CSR_PRECOND)\%crsptr\end{verbatim}
 if desired. To do that, use
  \begin{verbatim}FREEUPCSR(Elli\%G(CSR_PRECOND)\%crsptr)\end{verbatim}
 to free up the storage after the
  call to \begin{verbatim}CREATE_OPERATORS\end{verbatim}
 and recreate the preconditioner.
\item Iterate the solution. First, set the boundary conditions in
  Elli\%bc(i)\%p. Then call
  \begin{verbatim}take_div(ElliInfo,.true.).\end{verbatim}
 This routines
  calculates the divergence stored in div from the values in qe. The second
  arguments direct the routines to apply also the inhomogenous boundary
  conditions. This is done in two steps. First, the NM conditions on the
  velocity are applied directly on qe (i.e., on the edges the velocities are
  replaced with the desired ones). Then, along DR boundaries the gradient of a
  ``fake'' pressure obtained setting the ghost points of an otherwise zero
  pressure to their value. On the resulting field the divergence is taken.

Next, call
 \begin{verbatim}err=POISSON_SOLVER(Elli,SOLVE,solver),\end{verbatim}
 which will actually take
the step of solving the problem. 

Finally, call project(Elli,.true.) to apply the gradient of the potential
to the fields in eq. The second argument to true directs gradient to
add to the ghost point the component that comes from the interior (i.e., the
homogenous part).  At this stage, the values of the edge velocities on the NM
boundaries are not touched. If the metric contains off-diagonal terms, in view
of the aproximate nature of how NM boundary conditions are applied, there will
be an error in the divergence of the projected field. Thus one or more
iterations of the present step may be in order to reduce it to machine
precision. See how it is done in the examples.

\item When everything is done, and no more projections steps are expected,
  then the user calls
  \begin{verbatim}POISSON_SOLVERS(Elli,FINISH,solver)\end{verbatim}
 to release the
  memory associated to the LU decompositions, and a call to
  \begin{verbatim}FREEUP_OPERATORS(Elli)\end{verbatim}
 to release the memory taken up by the
  operators. 
\end{enumerate}

\section{Metric coefficients}

Up to this point, we have not discussed how the actual metric coefficients are
passed to the kit. It is up to the user to define then in a module that must
be called metric.f90 and must adhere to the format given in the mock-up module
in src. The kit needs a function g(i,j,k,l,p) which returns the metric
coefficient $g^{l,p}$ at location i,j,k. This must be defined as an elemental
function. When $l\neq p$ g refers to the metric on cell centers, otherwise is
the metric on cell edges on the appropriate edge. 
Also a logical function \begin{verbatim}g_mask(l,p)\end{verbatim} 
must be defined set to true if
$g^{l,p}\neq0$, false otherwise.     Again, it is imperative that the format
of g and \begin{verbatim}g_mask\end{verbatim} 
be as in the mockup module. Of course, other routines can be
put in the module, see examples. 

\section{2D vs. 3D}

In 2D, set mX(3)=1. 

\section{quirks}
\begin{itemize}
\item Because NM boundary conditions are implemented on the velocity, the potential
is not defined on the ghost cells on this boundary. 
\item qe actually stores the contravariant components multiplied by the
  jacobian of the transformation $\sqrt{g}$. Likewise, g is $\sqrt{g}\,g^{ij}$.
\end{itemize}

\end{document}


 


      


   
        

